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Abstract 



This is the manuscript of the chapter for a planned Handbook of Mathematical 
Methods in Imaging that surveys the mathematical models, problems, and algorithms 
of the Thermoacoustic (TAT) and Photoacoustic (PAT) Tomography. TAT and PAT 
represent probably the most developed of the several novel "hybrid" methods of medical 
imaging. These new modalities combine different physical types of waves (electromag- 
netic and acoustic in case of TAT and PAT) in such a way that the resolution and 
contrast of the resulting method are much higher than those achievable using only 
acoustic or electromagnetic measurements. 
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1 Introduction 



We provide here just a very brief description of the TAT/PAT procedure, since the relevant 
physics and biology details can be found in another chapter [132] in this volume, as well as 
in the surveys and books [133, 134, 137]. In TAT (PAT), a short pulse of radio- frequency EM 
wave (correspondingly, laser beam) irradiates a biological object (e.g., in the most common 
application, human breast), thus causing small levels of heating. The resulting thermoelastic 
expansion generates a pressure wave that starts propagating through the object. The ab- 
sorbed EM energy and the initial pressure it creates are much higher in the cancerous cells 
than in healthy tissues (see the discussion of this effect in [132-134, 137]). Thus, if one could 
reconstruct the initial pressure f(x), the resulting TAT tomogram would contain highly use- 
ful diagnostic information. The data for such a reconstruction are obtained by measuring 
time-dependent pressure p(x, t) using acoustic transducers located on a surface S (we will 
call it the observation or acquisition surface) completely or partially surrounding the body 
(see Fig. 1). Thus, although the initial irradiation is electro-magnetic, the actual recon- 
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Figure 1: TAT/PAT procedure with a partially surrounding acquisition surface. 



struction is based on acoustic measurements. As a result, the high contrast is produced due 
to a much higher absorption of EM energy by cancerous cells (ultrasound alone would not 
produce good contrast in this case), while the good (sub-millimeter) resolution is achieved 
by using ultrasound measurements (the radio frequency EM waves are too long for high- 
resolution imaging). Thus, TAT, by using two types of waves, combines their advantages, 
while eliminating their individual deficiencies. 

The physical principle upon which TAT/PAT is based was discovered by Alexander Gra- 
ham Bell in 1880 [21] and its application for imaging of biological tissues was suggested a 
century later [23]. It began to be developed as a viable medical imaging technique in the 
middle of 1990s [75,98]. 

Some of the mathematical foundations of this imaging modality were originally developed 
starting in the 1990s for the purposes of the approximation theory [84, 85] (see [7, 78] for 
extensive reviews of the resulting developments), integral geometry ([48, Chapter 5], [50]), 
and sonar and radar [27,86,93]. 

One can find recent reviews of the physics, biology, and mathematics issues of TAT/PAT 
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in [4, 44, 45, 69, 78, 99, 100, 102, 114, 128, 131, 133, 134, 137]. 

TAT/PAT is just one, probably the most advanced at the moment, example of the several 
recently introduced hybrid imaging methods, which combine different types of radiation to 
yield high quality of imaging unobtainable by single-radiation modalities (e.g., see [11,12, 
46, 79, 134] for other examples). 

2 Mathematical models of TAT 

In this section, we describe the commonly accepted mathematical model of the TAT pro- 
cedure and the main mathematical problems that need to be addressed. Since for all our 
purposes PAT results in the same mathematical model (although the biological features that 
TAT and PAT detect are different; see details in the chapter [15]), we will refer to TAT only. 

2.1 Point detectors and the wave equation model 

We will mainly assume that point-like omni-directional ultrasound transducers, located 
throughout an observation (acquisition) surface S, are used to detect the values of the pres- 
sure p(y, t), where y G S is a detector location and t > is the time of the observation. We 
also denote by c{x) the speed of sound at a location x. Then, it has been argued, that the 
following model describes correctly the propagating pressure wave p(x, t) generated during 
the TAT procedure (e.g., [15, 33, 127, 132, 135]): 

p tt = c 2 (x)A x p, i>0,iGK 3 
P(x, 0) = f(x),p t (x,0) = 0. 

Here f(x) is the initial value of the acoustic pressure, which one needs to find in order to 
create the TAT image. In the case of a closed acquisition surface S, we will denote by Q the 
interior domain it bounds. Notice that in TAT the function f(x) is naturally supported inside 
Q. We will see that this assumption about the support of / sometimes becomes crucial for 
the feasibility of reconstruction, although some issues can be resolved even if / has non-zero 
parts outside the acquisition surface. 

The data obtained by the point detectors located on a surface S are represented by the 
function 

g(y,t):=p(y,t)foiyeS,t>0. (2) 

Fig. 2 illustrates the space-time geometry of (1). 

We will incorporate the measured data g into the system (1), rewriting it as follows: 

Ptt = c 2 (x)A x p, i>0,xel 3 

p(x,0) = /(x) jft (x,0) = (3) 

p\s = g{y,t), ( y ,t)eSxR + . 

Thus, the goal in TAT/PAT is to find, using the data g(y,t) measured by transducers, 
the initial value f(x) at t — of the solution p(x,t) of (3). 
We will use the following notation: 
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Figure 2: The observation surface S and the domain f2 containing the object to be imaged. 



Definition 1. We will denote by W the forward operator 



W:f(x)^g(y,t), 



(4) 



where f and g are described in (3). 
Remark 2. 

• The reader should notice that if a different type of detectors is used, the system (1) stays 
intact, while the measured data will be represented differently from (2) (see Section 2.4)- 
This will correspondingly influence the reconstruction procedures. 

• We can consider the same problem in the space M. n of any dimension, not just in 3D. 
This is not merely a mathematical abstraction. Indeed, in the case of the so called 
integrating line detectors (Section 2.4), one deals with the 2D situation. 

2.2 Acoustically homogeneous media and spherical means 

If the medium being imaged is acoustically homogeneous (i.e., c(x) equals to a constant, 
which we will assume to be equal to 1 in appropriate units), as it is approximately the case 
in breast imaging, one deals with the constant coefficient wave equation problem 



In this case, the well known Poisson-Kirchhoff formulas [30, Ch. VI, Section 13.2, Formula 
(15)] for the solution of the wave equation gives in 3D: 





p(x,t) 



a-(t(Rf)(x,t)) , 



(6) 



where 




(7) 



is the spherical mean operator applied to the function f(x), dA is the standard area element 
on the unit sphere in IR 3 , and a is a constant. (Versions in all dimensions are known, see (16) 
and (15).) One can derive from here that knowledge of the function g(x, t) for x G S and all 
t > is equivalent to knowing the spherical mean Rf(x,t) of the function / for any points 
x G S and any t > 0. One thus needs to study the spherical mean operator R : f —+ Rf, or, 
more precisely, its restriction to the points x G S only, which we will denote by Ai: 



Due to the connection between the spherical mean operator and the wave equation, one can 
choose to work with the former, and in fact many works on TAT do so. The spherical mean 
operator Ai resembles the classical Radon transform, the common tool of computed tomog- 
raphy [73,88,89], which integrates functions over planes rather than spheres. This analogy 
with Radon transform, although often purely ideological, rather than technical, provides im- 
portant intuition and frequently points in reasonable directions of study. However, when the 
medium cannot be assumed to be acoustically homogeneous, and thus c(x) is not constant, 
the relation between TAT and integral geometric transforms, such as Radon transform or 
spherical mean, to a large extent breaks down, and thus one has to work with the wave 
equation directly. 

In what follows, we will address both models of TAT (the PDE model and the integral 
geometry model) and thus will deal with both forward operators W and Ai. 

2.3 Main mathematical problems arising in TAT 

We now formulate a list of problems related to TAT which will be addressed in detail in the 
rest of the article. (This list is more or less standard for a tomographic imaging method.) 

Sufficiency of the data. The first natural question to ask is: Is the data collected on 
the observation surface S sufficient for the unique reconstruction of the initial pressure 
f(x) (3)? In other words, is the kernel of the forward operator W zero? Or, to 
put it differently, for which sets S G M 3 the data collected by transducers placed 
along S determines / uniquely? Yet another interpretation of this question is through 
observability of solutions of the wave equation on the set S: does observation on S of 
a solution of the problem (1) determine the solution uniquely? 

When the speed of sound is constant, and thus the spherical mean model applies, the 
equivalent question is whether the operator Ai has zero kernel on an appropriate class 
of functions (say, continuous functions with compact support) 

As it is explained in [7], the choice of precise conditions on the local function class, 
such as continuity, is of no importance for the answer to the uniqueness question, while 
behavior at infinity (e.g., compactness of support) is. So, without loss of generality, 
when discussing uniqueness, one can assume f(x) in (3) to be infinitely differentiable. 

Inversion formulas and algorithms. Since a practitioner needs to see the actual tomo- 
gram, rather than just know its existence, the next natural question arises: If unique- 
ness the data collected on S is established, what are the actual inversion formulas or 
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algorithms? Here again one can work with smooth functions, in the end extending the 
formulas by continuity to a wider class. 

Stability of reconstruction. If we can invert the transform and reconstruct / from the 
data g, how stable is the inversion? The measured data are unavoidably corrupted by 
errors, and stability means that small errors in the data lead to only small errors in 
the reconstructed tomogram. 

Incomplete data problems. What happens if the data is "incomplete," for instance 
if one can only partially surround the object by transducers? Does this lead to any 
specific deterioration in the tomogram, and if yes, to what kind of deterioration? 

Range descriptions. 

The next question is known to be important in analysis of tomographic problems: 
What is the range of the forward operator W : / i— > g that maps the unknown function 
/ to the measured data gl In other words, what is the space of all possible "ideal" data 
g(t, y) collected on the surface SI In the constant speed of sound case, this is equivalent 
to the question of describing the range of the spherical mean operator M. in appropriate 
function spaces. Such ranges often have infinite co-dimensions, and the importance of 
knowing the range of Radon type transforms for analyzing problems of tomography is 
well known. For instance, such information is used to improve inversion algorithms, 
complete incomplete data, discover and compensate for certain data errors, etc. (e.g., 
[35,47-49,63-65,76,88,89,102] and references therein). In TAT, range descriptions 
are also closely connected with the speed of sound determination problem listed next 
(see Section 3.6 for a discussion of this connection). 

Speed of sound reconstruction. As the reader can expect, reconstruction procedures 
require the knowledge of the speed of sound c(x). Thus, the problem arises of the 
recovery of c(x) either from an additional scan, or (preferably) from the same TAT 
data. 

2.4 Variations on the theme: planar, linear, and circular integrat- 
ing detectors 

In the described above most basic and well-studied version of TAT, one utilizes point-like 
broadband transducers to measure the acoustic wave on a surface surrounding the object 
of interest. The corresponding mathematical model is described by the system (3). In 
practice, the transducers cannot be made small enough, since smaller detectors yield weaker 
signals resulting in low signal-to-noise ratios. Smaller transducers are also more difficult to 
manufacture. 

Since finite size of the transducers limits the resolution of the reconstructed images, 
researchers have been trying to design alternative acquisition schemes using receivers that 
are very thin but long or wide. Such are 2D planar detectors [25,56] and ID linear and 
circular [26, 52, 106, 142] detectors. 

We will assume throughout this section that the speed of sound c(x) is constant and 
equal to 1. 
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Planar detectors are made from a thin piezoelectric polymer film glued onto a flat sub- 
strate (see, for example [109]). Let us assume that the object is contained within the sphere 
of radius R. If the diameter of the planar detector is sufficiently large (see [109] for details), 
it can be assumed to be infinite. The mathematical model of such an acquisition technique 
is no longer described by (3). Let us define the detector plane H(s, uj) by equation x ■ uj — s, 
where u is the unit normal to the plane and s is the (signed) distance from the origin to the 
plane. Then, while the propagation of acoustic waves is still modeled by (1), the measured 
data g p ianar(s, t, uj) (up to a constant factor which we will, for simplicity, assume to be equal 
to 1) can be represented by the following integral: 



g P lanar(s,UJ,t) = j p(x,t)dA(x) 

n(s,w) 

where dA(x) is the surface measure on the plane. Obviously, 

9planar(s,UJ,0) = J p(x,0)dA(x) = J f (x)dA(x) = F (s, Uj) , 

n(s,w) n(s,£j) 

i.e. the value of g at t = coincides with the integral F(s,uj) of the initial pressure f(x) 
over the plane H(s,uj) orthogonal to uj. 

One can show [25,56] that for a fixed uj, function g p i anar {s,u,t) is the solution to ID 
wave equation 

d 2 g d 2 g 



ds 2 dt 2 ' 



and thus 



gplanar(s, UJ , t) — — [g p lanar{s, UJ , S — t) + g p l a nar(s, UJ , S + t)] 

= -[F(s + t,w)+F(s-t,u)]. 

Since the detector can only be placed outside the object, i.e. s > R, the term F(s + t,u) 
vanishes, and one obtains 

g P lanar{s, UJ , t) = F(s -t,u). 

In other words, by measuring g p i an ar(s, uj, t), one can obtain values of the planar integrals of 
f(x). If, as proposed in [25,56], one conducts measurements for all planes tangent to the 
upper half-sphere of radius R (i.e. s = R, uj 6 S+), then the resulting data yield all values 
of the standard Radon transform of f{x). Now the reconstruction can be carried out using 
one of the many known inversion algorithms for the latter transform (see [73, 88, 89]). 

Linear detectors are based on optical detection of acoustic signal. Some of the proposed 
optical detection schemes utilize as the sensitive element a thin straight optical fiber in com- 
bination with Fabry- Perot interferometer [26,52]. Changes of acoustic pressure on the fiber 
change (proportionally) its length; this elongation, in turn, is detected by interferometer. A 
similar idea is used in [106]; in this work the role of a sensitive element is played by a laser 
beam passing through the water in which the object of interest is submerged, and thus the 
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measurement does not perturb the acoustic wave. In both cases, the length of the sensitive 
element exceeds the size of the object, while the diameter of the fiber (or of the laser beam) 
can be made extremely small (see [109] for a detailed discussion), which removes restrictions 
on resolution one can achieve in the images. 

Let us assume that the fiber (or laser beam) is aligned along the line l(s\, s 2 , u>i, 0J2) = 
{x\x = S\Ui + 52^2 + su}, where vectors u)i,u 2 , and uj form an ortho-normal basis in R 3 . 
Then the measured quantities gu n ear(si, S2, 0J\, 0J2, t) are equal (up to a constant factor which, 
we will assume, equals to 1) to the following line integral: 



Similarly to the case of planar detection, one can show [26,52,106], that for fixed vectors 
loi,uj2 the measurements gu near (si , s 2 , uj 2 , t) satisfy the 2D wave equation 



The initial values gu n ear(si, s 2 , c^i, OJ21 0) coincide with the line integrals of f(x) along lines 
l(si, S2, U)i, UJ2) ■ Suppose one makes measurements for all values of s\(t), s 2 (t) corresponding 
to a curve 7 = {x\x = Si(t)u)i + s 2 (r)u;2, t~o < t < r\} lying in the plane spanned by loi, uj 2 . 
Then one can try to reconstruct the initial value of g from the values of g on 7. This problem 
is a 2D version of (3) and thus the known algorithms (see Section 4) are applicable. 

In order to complete the reconstruction from data obtained using line detectors, the 
measurements should be repeated with different directions of u. For each value of u the 2D 
problem is solved; the solutions of these problems yield values of line integrals of f(x). If 
this is done for all values of u lying on a half circle, the set of the recovered line integrals 
of f(x) is sufficient for reconstructing this function. Such a reconstruction represents the 
inversion of the well known in tomography X-ray transform. The corresponding theory and 
algorithms can be found, for instance, in [73,88,89]. 

Finally, the use of circular integrating detectors was considered in [142]. Such a detector 
can be made out of optical fiber combined with an interferometer. In [142], a closed form 
solution of the corresponding inverse problem is found. However, this approach is very new 
and neither numerical examples, nor reconstructions from real data have been obtained yet. 

3 Mathematical analysis of the problem 

In this section, we will address most of the issues described in Section 2.3, except the recon- 
struction algorithms, which will be discussed in Section 4. 

3.1 Uniqueness of reconstruction 

The problem discussed here is the most basic one for tomography: given an acquisition 
surface S along which we distribute detectors, is the data g(y,t) for y e S, t > (see (3)) 
sufficient for a unique reconstruction of the tomogram /? A simple counting of variables 




d 2 g d 2 g d 2 g 
ds\ ds\ dt 2 
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shows that S should be a hyper-surface in the ambient space (i.e., a surface in M 3 or a 
curve in M 2 ). As we will see below, although there are some simple counter-examples and 
remaining open problems, for all practical purposes, the uniqueness problem is positively 
resolved, and most surfaces S do provide uniqueness. We address this issue for acoustically 
homogeneous media first and then switch to the variable speed case. 

Before doing so, however, we would like to dispel a concern that arises when one looks 
at the problem of recovering / from g in (3). Namely, an impression might emerge that 
we consider an initial-boundary value (IBV) problem for the wave equation in the cylinder 
Q x M + , and the goal is to recover the initial data / from the known boundary data g. This 
is clearly impossible, since according to standard PDE theorems (e.g., [30, 36]), one can solve 
this IBV problem for arbitrary choice of the initial data / and boundary data g (as long 
as they satisfy simple compatibility conditions, which are fulfilled for instance if / vanishes 
near S and g vanishes for small t, which is the case in TAT). This means that apparently g 
contains essentially no information about / at all. This argument, however, is flawed, since 
the wave equation in (3) holds in the whole space, not just in Q. In other words, S is not a 
boundary, but rather an observation surface. In particular, considering the wave equation in 
the exterior of S, one can derive that if / is supported inside Q, the boundary values g of the 
solution p of (3) also determine the normal derivative of p at S for all positive times. Thus, 
we in fact have (at least theoretically) the full Cauchy data of the solution p on S, which 
should be sufficient for reconstruction. Another way of addressing this issue is to notice that 
if the speed of sound is constant, or at least non-trapping (see the definition below in Section 
3.1.2), the energy of the solution in any bounded domain (in particular, in Q) must decay in 
time. The decay when t — > oo together with the boundary data g guarantee the uniqueness 
of solution, and thus uniqueness of recovery /. 

These arguments, as the reader will see, play a role in understanding reconstruction 
procedures. 

3.1.1 Acoustically homogeneous media 

We assume here the sound speed c(x) to be constant (in appropriate units, one can choose 
it to be equal to 1, which we will do to simplify considerations). 

In order to state the first important result on uniqueness, let us recall the system (5), 
allowing an arbitrary dimension n of the space: 



We introduce the following useful definition: 

Definition 3. A set S is said to be uniqueness set, if when used as the acquisition surface, 
it provides sufficient data for unique reconstruction of the compactly supported tomogram f 
(i.e., the observed data g in (9) determines uniquely function f). Otherwise, it is called a 
non-uniqueness set. 

In other words, S is a uniqueness set if the forward operator W (or, equivalently, M) has 
zero kernel. 




(9) 



10 



We will start with a very general statement about the acquisition (observation) sets S 
that provide insufficient information for unique reconstruction of / (see [7] for the proof and 
references) : 

Theorem 4. If S is a non-uniqueness set, then there exists a non-zero harmonic polynomial 
Q, which vanishes on S. 

This theorem implies, in particular, that all "bad" (non-uniqueness) observation sets 
are algebraic, i.e. have a polynomial vanishing on them. Turning this statement around, 
we conclude that any set S that is a uniqueness set for harmonic polynomials, is sufficient 
for unique TAT reconstruction (although, as we will see in Section 3.3, this does not mean 
practicality of the reconstruction). 

The proof of Theorem 4, which the reader can find in [7, 78], is not hard and in fact is 
enlightening, but providing it would lead us too far from the topic of this survey. 

We will consider first the case of closed acquisition surfaces, i.e. the ones that completely 
surround the object to be imaged. We will address the general situation afterwards. 

Closed acquisition surfaces S 

Theorem 5. ([7]) If the acquisition surface S is the boundary of bounded domain Q (i.e., 
a closed surface), then it is a uniqueness set. Thus, the observed data g in (9) determines 
uniquely the sought function f G L^ omp (M. n ). (The statement holds, even though f is not 
required to be supported inside S .) 

Proof: Indeed, since there are no non-zero harmonic functions vanishing on a closed 
surface S, Theorem 4 implies Theorem 5. □ 

There is, however, another, more intuitive, explanation of why Theorem 5 holds true 
(although it requires somewhat stronger assumptions, or a more delicate proof than the one 
indicated below). Namely, since the solution p of (9) has compactly supported initial data, 
its energy is decaying inside any bounded domain, in particular inside Q (see Section 3.1.2 
and [34,67] and references therein about local energy decay). On the other hand, if there 
is non-uniqueness, there exists a non-zero / such that g(y,t) = for all y G S and t. This 
means that we can add homogeneous Dirichlet boundary conditions p \s— to (9). But then 
the standard PDE theorems [30, 36] imply that the energy stays constant in Q. Combination 
of the two conclusions means that p is zero in Q for all times t. It is well known [30] that 
such a solution of the wave equation must be identically zero everywhere, and thus / = 0. 

This energy decay consideration can be extended to some classes of non-compactly sup- 
ported functions / of the IP classes, leading to the following result of [1]: 

Theorem 6. [1] Let S be the boundary of a bounded domain in M n and / G L p (M n ). Then 

1- Ifp < ~[ an d the spherical mean of f over almost every sphere centered on S is equal 
to zero, then f = 0. 

2. The previous statement fails when p > and S is a sphere. 

In other words, a closed surface S is a uniqueness set for functions f G L p (W l ) whenp < — 
and might fail to be such when p > . 



11 



This result shows that the assumption, if not necessarily of compactness of support of 
/, but at least of a sufficiently fast decay of / at infinity, is important for the uniqueness to 
hold. 

General acquisition sets S 

Theorems 4 and 5 imply the following useful statement: 

Theorem 7. If a set S is not algebraic, or if it contains an open part of a closed analytic 
surface T, then it is a uniqueness set. 

Indeed, the first claim follows immediately from Theorem 4. The second one works out 
as follows: if an open subset of an analytic surface T is a non-uniqueness set, then by an 
analytic continuation type argument (see [7]), one can show that the whole T is such a set. 
However, this is impossible, due to Theorem 5. 

There are simple examples of non-uniqueness surfaces. Indeed, if S is a plane in 3D 
(or a line in 2D, or a hyperplane in dimension n) and f(x) in (3) is odd with respect to 
S, then clearly the whole solution of (3) has the same parity and thus vanishes on 5* for 
all times t. This means that, if one places transducers on a planar S, they might register 
zero signals at all times, while the function / to be reconstructed is not zero. Thus, there 
is no uniqueness of reconstruction when S is a plane. On the other hand (see [30,72]), if 
/ is supported completely on one side of the plane S (the standard situation in TAT), it 
is uniquely recoverable from its spherical means centered on S, and thus from the observed 
data g. 

The question arises what are other "bad" (non-uniqueness) acquisition surfaces than 
planes. This issue has been resolved in 2D only. Namely, consider a set of N lines on the 
plane intersecting at a point and forming at this point equal angles. We will call such a 
figure the Coxeter cross Etv (see Fig. 3). it is easy to construct a compactly supported 




Figure 3: Coxeter cross of N lines. 

function that is odd simultaneously with respect of all lines in E^. Thus, a Coxeter cross is 
also a non-uniqueness set. The following result, conjectured in [84, 85] and proven in the full 
generality in [7], shows that, up to adding finitely many points, this is all that can happen 
to non-uniqueness sets: 
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Theorem 8. [7] A set S in the plane M 2 is a non-uniqueness set for compactly supported 
functions f , if and only if it belongs to the union Sjy (J $ of a Coxeter cross Sjy and a finite 
set of points $. 

Again, compactness of support is crucial for the proof provided in [7] . There are no other 
proofs known at the moment of this result (see the corresponding open problem in Section 
5). In particular, there is no proven analog of Theorem 6 for non-closed sets S (unless S is 
an open part of a closed analytic surface). 

The n-dimensional (in particular, 3D) analog of Theorem 8 has been conjectured [7], but 
never proven, although some partial advances in this direction have been made in [8,42]. 

Conjecture 9. A set S in IR n is a non-uniqueness set for compactly supported functions f , 
if and only if it belongs to the union S [J $, where £ is the cone of zeros of a homogeneous 
(with respect to some point in M. n ) harmonic polynomial, and $ is an algebraic sub-set ofW a 
of dimension at most n — 2 (see Fig. 4)- 




Figure 4: The conjectured structure of a most general non-uniqueness set in 3D. 

Uniqueness results for a finite observation time 

So far, we have addressed only the question of uniqueness of reconstruction in the non- 
practical case of the infinite observation time. There are, however, results that guarantee 
uniqueness of reconstruction for a finite time of observation. The general idea is that it is 
sufficient to observe for the time that it takes the geometric rays (see Section 3.1.2) from 
the interior Q of if? to reach S. In the case of a constant speed, which we will assume to be 
equal to 1, the rays are straight and are traversed with the unit speed. This means that if 
D is the diameter of Q (i.e., the maximal distance between two points in the closure of Q), 
then after time t = D, all rays coming from Q have left the domain. Thus, one hopes that 
waiting till time t = D might be sufficient. In fact, due to the specific initial conditions in 
(3), namely, that the time derivative of the pressure is equal to zero at the initial moment, 
each singularity of / emanates two rays, and at least one of them will reach S in time not 
exceeding D/2. And indeed, the following result of [42] holds: 

Theorem 10. [42] If S is smooth and closed surface bounding domain Q and D is the 
diameter of Q, then the TAT data on S collected for the time < t < 0.5D, uniquely 
determines f . 
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Notice that a shorter collection time does not guarantee uniqueness. Indeed, if S is a 
sphere and the observation time is less than 0.5D, due to the finite speed of propagation, no 
information from a neighborhood of the center can reach S during observation. Thus, values 
of / in this neighborhood cannot be reconstructed. 

3.1.2 Acoustically inhomogeneous media 

We assume that the speed of sound is strictly positive, c(x) > c > 0, and such that c(x) — 1 
has compact support, i.e. c(x) = 1 for large x. 

Trapping and non-trapping 

We will frequently impose the so-called non-trapping condition on the speed of sound c(x) 
in R n . To introduce it, let us consider the Hamiltonian system in R^l with the Hamiltonian 



The solutions of this system are called bicharacteristics and their projections into M" are 
rays (or geometric rays). 

Definition 11. We say that the speed of sound c(x) satisfies the non-trapping condition, if 
all rays with £o 7^ tend to infinity when t — > 00. 

The rays that do not tend to infinity, are called trapped. 

A simple example, where quite a few rays are trapped, is the radial parabolic sound speed 
c(x) = c\x\ 2 . 

It is well known (e.g., [66]) that singularities of solutions of the wave equation are carried 
by geometric rays. In order to make this statement more precise, we need to recall the notion 
of a wave front set WF(u) of a distribution u(x) in K n . This set carries detailed information 
on singularities of u(x). 

Definition 12. Distribution u(x) is said to be microlocally smooth near a point (xo>£o)> 
where xq,C,o £ E n and^Q 7^ 0, if there is a smooth 1 'cut- off" function 4>(x) such that 4>(xq) 7^ 
and that the Fourier transform </>«(£) of the function (p(x)u(x) decays faster than any power 
\£,\~ N when |£| — > 00, in directions that are close to the direction of £o- 1 

The wave front set WF(u) C R™ x (R£ \ 0) of u consists of all pairs (x , £0) suc h that 
u is not microlocally smooth near (xo,£o)- 

In other words, if (xo,£o) WF(u), then u is not smooth near Xq, and the direction of 
£0 indicates why it is not: the Fourier transform does not decay well in this direction. For 
instance, if u{x) consists of two smooth pieces joined non-smoothly across a smooth interface 
S, then WF(u) can only contain pairs (x, £) such that 16E and £ is normal to S at x. 

1 We remind the reader that if this Fourier transform decays that way in all directions, then u{x) is 
smooth (infinitely differcntiable) near the point xq. 



H= c ^f±\Z\ 2 : 




(10) 
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It is known that the wave front sets of solutions of the wave equation propagate with 
time along the bicharacteristics introduced above. This is a particular instance of a more 
general fact that applies to general PDEs and can be found in [66, 121]. As a result, if after 
time T all the rays leave the domain Q of interest, the solution becomes smooth (infinitely 
differentiable) inside Q. 

One can find simple introduction to the notions of microlocal analysis, such as the wave 
front set, for instance in [126], and more advanced versions in [66,121]. Applications of 
microlocal analysis to integral geometry are discussed in [51,53-55]. 

The notion of so called local energy decay, which we survey next, is important for the 
understanding of the non-trapping conditions in TAT. 

Local energy decay estimates 

Assuming that the initial data f(x) (1) is compactly supported and the speed c(x) is 
non-trapping, one can provide the so called local energy decay estimates [34,129,130]. 
Namely, in any bounded domain Q, the solution p(x, t) of (1) satisfies, for a sufficiently large 
To and for any (k,m), the estimate 

< C kjm u k {t)\\f\\ L 2, for i£fl,t> T . (11) 

Here Vk(t) = t~ n+1 ~ k for even n and v k {t) = e~ 5t for odd n and some 5 > 0. Any value To 
larger than the diameter of Q works in this estimate. 

Uniqueness result for non-trapping speeds 

If the speed is non-trapping, the local energy decay allows one to start solving the problem 
(3) from t = oo, imposing zero conditions at t = oo and using the measured data g as the 
boundary conditions. This leads to recovery of the whole solution, and in particular its 
initial value f(x). As the result, one obtains the following simple uniqueness result of [3]: 

Theorem 13. [3] If the speed c(x) is smooth and non-trapping and the acquisition surface 
S is closed, then the TAT data g(y,t) determines the tomogram f(x) uniquely. 

Notice that the statement of the theorem holds even if the support of / is not completely 
inside of the acquisition surface S. 

Uniqueness results for finite observation times 

As in the case of constant coefficients, if the speed of sound is non-trapping, appropriately 
long finite observation time suffices for the uniqueness. Let us denote by T(Q) the supremum 
of the time it takes the ray to reach S, over all rays originating in Q. In particular, if c(x) 
is trapping, T(Q) might be infinite. 

Theorem 14. [123] The data g measured till any time T larger than T(Q) is sufficient for 
unique recovery of f . 



Qk+\m 

d k t df 
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3.2 Stability 



By stability of reconstruction of the TAT tomogram / from the measured data g we mean 
that small variations of g in an appropriate norm lead to small variations of the reconstructed 
tomogram /, also measured by an appropriate norm. In other words, small errors in the 
data lead to small errors in the reconstruction. 

We will try to give the reader a feeling of the general state of affairs with stability, 
referring to the literature (e.g., [5,68, 78, 103, 123]) for further exact details. 

We will consider as functional spaces the standard Sobolev spaces H s of smoothness s. 
We will also denote, as before, by W the operator transforming the unknown / into the data 
9- 

Let us recall the notions of Lipschitz and Holder stability. An even weaker loga- 
rithmic stability will not be addressed here. The reader can find discussion of the general 
stability notions and issues, as applied to inverse problems, in [70]. 

Definition 15. The operation of reconstructing f from g is said to be Lipschitz stable 

between the spaces H S2 and H Sl , if the following estimate holds for some constant C: 

< C\\g\\ H s 2 . 

The reconstruction is said to be Holder stable (a weaker concept), if there are constants 
s i) s 2, S3, C, ji > 0, and 5 > such that 

\\f\\Hn<C\\gf H32 

for all f such that \\f\\H s 3 < te- 
stability can be also interpreted in the terms of the singular values aj of the forward 
operator / 1— > g in L 2 , which have at most power decay when j — > 00. The faster is the 
decay, the more unstable the reconstruction becomes. The problems with singular values 
decaying faster than any power of j are considered to be extremely unstable. Even worse 
are the problems with exponential decay of singular values (analytic continuation or solving 
Cauchy problem for an elliptic operator belong to this class). Again, the book [70] is a good 
source for finding detailed discussion of such issues. 

Consider as an example inversion of the standard in X-ray CT and MRI Radon transform 
that integrates a function / over hyper-planes in M n . It smoothes function by "adding 
(n — l)/2 derivatives." Namely, it maps continuously if s -functions in Q into the Radon 
projections of class _£/" s +( n_1 )/ 2 . Moreover, the reconstruction procedure is Lipshitz stable 
between these spaces (see [88] for detailed discussion). 

One should notice that since the forward mapping is smoothing (it "adds derivatives" 
to a function), the inversion should produce functions that are less smooth than the data, 
which is an unstable operation. The rule of thumb is that the stronger is smoothing, the 
less stable is inversion (this can be rigorously recast in the language of the decay of singular 
values). Thus, problems that require reconstructing non-smooth functions from infinitely 
different iable (or even worse, analytic) data, are extremely unstable (with super-algebraic 
or exponential decay of singular values correspondingly). This is just a consequence of the 
standard Sobolev embedding theorems (see, e.g., how this applies in TAT case in [91]). 
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In the case of a constant sound speed and the acquisition surface completely surrounding 
the object, as we have mentioned before, the TAT problem can be recast as inversion of the 
spherical mean transform M. (see Section 2). Due to analogy between the spheres centered 
on S and hyperplanes, one suspects that inversion of the spherical mean operator 
Ai is as Lipschitz stable as the inversion of the Radon transform. This indeed is 
the C3iSG ; clS long as / is supported inside S, as can be found in [103]. In the cases 
when closed form inversion formulas are available (see Section 4.1.1), this stability can also 
be extracted from them. If the support of / does reach outside, reconstruction of the 
part of / that is outside is unstable (i.e., is not even Holder stable, due to the reasons 
explained in Section 3.3). 

In the case of variable non-trapping speed of sound c(x), integral geometry does 
not apply anymore, and one needs to address the issue using, for instance, time reversal. In 
this case, stability follows by solving the wave equation in reverse time starting from t — oo, 
as it is done in [3]. In fact, Lipschitz stability in this case holds for any observation 
time exceeding T(Q) (see [123], where microlocal analysis is used to prove this result). 

The bottom line is that TAT reconstruction is sufficiently stable, as long as the 
speed of sound is non-trapping. 

However, trapping speed does cause instability [68]. Indeed, since some of the rays are 
trapped inside Q, the information about some singularities never reaches S (no matter for 
how long one collects the data), and thus, as it is shown in [91], the reconstruction is not 
even Holder stable, and the singular values have super-algebraic decay. See also Section 3.3 
below for a related discussion. 

3.3 Incomplete data 

In the standard X-ray CT, incompleteness of data arises, for instance, if not all projection 
angles are accessible, or irradiation of certain regions is avoided, or as in the ROI (region of 
interest) imaging, only the ROI is irradiated. 

It is not that clear what incomplete data means in TAT. Usually one says that one 
deals with incomplete TAT data, if the acquisition surface does not surround the 
object of imaging completely. For instance, in breast imaging it is common that only a 
half-sphere arrangement of transducers is possible. We will see, however, that incomplete 
data effects in TAT can also arise due to trapping, even if the acquisition surface 
completely surrounds the object. 

The questions addressed here are: 

1. Is the collected incomplete data sufficient for unique reconstruction? 

2. If yes, does the incompleteness of the data have any effect on stability and quality 
of the reconstruction? 

3.3.1 Uniqueness of reconstruction 

Uniqueness of reconstruction issues can be considered essentially resolved for incomplete data 
in TAT, at least in most situations of practical interest. We will briefly survey here some of 
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the available results. In what follows, the acquisition surface S is not closed (otherwise the 
problem is considered to have complete data). 

Uniqueness for acoustically homogeneous media 

In this case, Theorem 7 contains some useful sufficient conditions on S that guarantee 
uniqueness. Microlocal results of [7,86,122], as well as the PDE approach of [42] further 
applied in [8] provide also some other conditions. We assemble some of these in the following 
theorem: 

Theorem 16. Let S be a non-closed acquisition surface in TAT. Each of the following 
conditions on S is sufficient for the uniqueness of reconstruction of any compactly supported 
function f from the TAT data collected on S : 

1. Surface S is not algebraic (i.e., there is no non-zero polynomial vanishing on S). 

2. Surface S is a uniqueness set for harmonic polynomials (i.e., there is no non-zero 
harmonic polynomial vanishing on S). 

3. Surface S contains an open piece of a closed analytic surface T. 

4- Surface S contains an open piece of an analytic surface T separating the space M. n such 
that f is supported on one side ofT. 

5. For some point y £ S, the function f is supported on one side of the tangent plane T y 
to S at y. 

For instance, if the acquisition surface S is just a tiny non-algebraic piece of a surface, 
data collected on S determines the tomogram / uniquely. However, one realizes that such 
data is unlikely to be useful for any practical reconstruction. Here the issue of stability of 
reconstruction kicks in, as it will be discussed in the stability sub-section further down. 

Uniqueness for acoustically inhomogeneous media 

In the case of a variable speed of sound, there still are uniqueness theorems for partial 
data [123,124], e.g. 

Theorem 17. [123] Let S be an open part of the boundary dfl of a strictly convex domain 
Q and the smooth speed of sound equals 1 outside Q. Then the TAT data collected on S for a 
time T > T(Q) determines uniquely any function f £ Hq(Q), whose support does not reach 
the boundary. 

A modification of this result that does not require strict convexity is also available in 
[124]. 

While useful uniqueness of reconstruction results exist for incomplete data problems, all 
such problems are expected to show instability. This issue is discussed in the sub-sections 
below. This will also lead to a better understanding of incomplete data phenomena in TAT. 
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3.3.2 "Visible" ("audible") singularities 

According to the discussion in Section 3.1.2, the singularities (the points of the wave front 
set WF(f) of the function / in (3)) are transported with time along the bi-characteristics 
(10). Thus, in the x-space they are transported along the geometric rays. These rays may 
or may not reach the acquisition surface S, which triggers the introduction of the following 
notion: 

Definition 18. A phase space point (xo, £o) ^ s sa ^d to be "visible" (sometimes the word 
"audible" is used instead), if the corresponding ray (see (10)) reaches in finite time the 
observation surface S. 

A region U C W 1 is said to be in the visibility zone, if all points (xq, £o) with xq G U 
are visible. 

An example of wave propagation through inhomogeneous medium is presented in Figure 
5. The open observation surface S in this example consists of the two horizontal and the 
left vertical sides of the square. Figure 5(a) shows some rays that bend, due to acoustic 
inhomogeneity, and leave through the opening of the observation surface S (the right side 
of the square). Fig. 5 (b) presents a flat phantom, whose wavefront set creates these 
escaping rays, and thus is mostly invisible. Then Fig. 5 (c-f) show the propagation of the 
corresponding wave front. 




(d) (e) (f) 

Figure 5: (a) Some rays starting along the interval x G [—0.7, —0.2] in the vertical directions 
escape on the right; (b) a flat phantom with "invisible wavefront"; (c-f)propagation of the 
flat front: most of the energy of the signal leaves the square domain through the hole on the 
right. 

Since the information about the horizontal boundaries of the phantom escapes, one does 
not expect to reconstruct it well. Fig. 6 shows two phantoms and their reconstructions 
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from the partial data: (a-b) correspond to the vertical flat phantom, whose only invisible 
singularities are at its ends. One sees essentially good reconstruction, with a little bit of 
blurring at the endpoints. On the other hand, reconstruction of the horizontal phantom 
with almost the whole wave front set invisible, does not work. The next Fig. 7 shows a more 




(a) (b) 




Figure 6: Reconstruction with the same speed of sound: (a-b) phantom with strong vertical 
fronts and its reconstruction; (c-d) phantom with strong horizontal fronts and its reconstruc- 
tion. 

complex square phantom, whose singularities corresponding to the horizontal boundaries 
are invisible, while the vertical boundaries are fine. One sees clearly that the invisible parts 
have been blurred away. On the other hand, Fig. 11(a) in Section 4 shows that one can 




(a) (b) (c) 

Figure 7: Reconstruction with the same speed of sound: (a) phantom; (b) its reconstruction; 
(c) a magnified fragment of (b). 
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reconstruct an image without blurring and with correct values, if the image is located in the 
visibility region. The reconstructed image in this figure is practically indistinguishable from 
the phantom shown in Figure 10(a). 

Remark 19. If S is a closed surface and xq is a point outside of S, there is a vector £o 7^ 
such that (xo,£o) "invisible." Thus, the visibility zone does not reach outside the closed 
acquisition surface S. 

3.3.3 Stability of reconstruction for incomplete data problems 

In all examples above, uniqueness of reconstruction held, but the images were still blurred. 
The question arises whether the blurring of "invisible" parts is avoidable (after all, the 
uniqueness theorems seem to claim that "everything is visible"). The answer to this is, in 
particular, the following result of [91], which is an analog of similar statements in X-ray 
tomography: 

Theorem 20. [91] If there are invisible points (xq, £0) in Q x (R^ \ 0), then inversion of the 
forward operator W is not Holder stable in any Sobolev spaces. The singular values o~j ofW 
in L 2 decay super- algebraically. 

Thus, having invisible singularities makes the reconstruction severely ill-posed. In par- 
ticular, according to Remark 19, this theorem implies the following statement: 

Corollary 21. Reconstruction of the parts of f(x) supported outside the closed observation 
surface S is unstable. 

On the other hand, 

Theorem 22. [123] All visible singularities of f can be reconstructed with Lipschitz stability 
(in appropriate spaces). 

Such a reconstruction of visible singularities can be obtained in many ways, for instance 
just by replacing the missing data by zeros (with some smoothing along the junctions with 
the known data, in order to avoid artifact singularities). However, there is no hope for stable 
recovery of the correct values of f(x), if there are invisible singularities. 

3.4 Discussion of the visibility condition 

Visibility for acoustically homogeneous media 

In the constant speed case, the rays are straight, and thus the visibility condition has a 
simple test: 

Proposition 23. (e.g., [68, 140, 141]) If the speed is constant, a point x$ is in the visible 
region, if and only if any line passing through xq intersects at least once the acquisition 
surface S (and thus a detector location). 

Figure 8 illustrates this statement. It shows a square phantom and its reconstruction 
from complete data and from the data collected on the half-circle S surrounding the left half 
of object. The parts of the interfaces where the normal to the interface does not cross S are 
blurred. 
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(a) (b) (c) 

Figure 8: Reconstruction from incomplete data using closed form inversion formula in 2D; 
detectors are located on the left half circle of radius 1.05 (a) phantom (b) reconstruction 
from complete data (c) reconstruction from the incomplete data 

Visibility for acoustically inhomogeneous media 

When the speed of sound is variable, an analog of Proposition 23 holds, with lines replaced 
by rays. 

Proposition 24. (e.g., [68, 91, 123]) A point xq is in the visible region, if and only if for 
an U £o 7^ at least one of the two geometric rays starting at (xo,£o) an d at (xo, — £o) (see 
(10)) intersects the acquisition surface S (and thus a detector location). 

The reader can now see an important difference between the acoustically homogeneous 
and inhomogeneous media. Indeed, even if S surrounds the support of / completely, trapped 
rays will never find their way to S, which will lead, as we know by now, to instabilities and 
blurring of some interfaces. 

Thus, presence of rays trapped inside the acquisition surface creates effects of incomplete 
data type. This is exemplified in Fig. 9 with a square phantom and its reconstruction shown 
in the presence of a trapping (parabolic) speed. Notice that the square centered at the center 
of symmetry of the speed is reconstructed very well (see (d)), since none of the rays carrying 
its singularities is trapped. 

3.5 Range conditions 

In this section we address the problem of describing the ranges of the forward operators W 
(see (4)) and M. (see (8)), the latter in the case of an acoustically homogeneous medium 
(i.e., for c = const). The ranges of these operators, similarly to the range of the Radon and 
X-ray transforms (see [88,89]), are of infinite co-dimensions. This means that ideal data g 
from a suitable function space satisfy infinitely many mandatory identities. Knowing the 
range is useful for many theoretical and practical purposes in various types of tomography 
(reconstruction algorithms, error corrections, incomplete data completion, etc.), and thus 
this topic has attracted a lot of attention (e.g., [35, 47-49, 63, 64, 76, 80, 88, 89, 102, 119] and 
references therein). 

As we will see in the next section, range descriptions in TAT are also intimately related 
to recovery of the unknown speed of sound. 
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(a) 



(b) 



(c) 



(d) 



Figure 9: Reconstruction of a square phantom from full data in the presence of a trapping 
parabolic speed of sound (the speed is radial with respect to the center of the picture): (a) an 
off-center phantom; (b) its reconstruction; (c) a magnified fragment of (b); (d) reconstruction 
of a centered square phantom. 

We recall [47, 48, 63, 88] that for the standard Radon transform 



where / is assumed to be smooth and supported in the unit ball B — {x \ \x\ < 1}, the 
range conditions on g(s,u) are: 

1. smoothness and support: g £ ([— 1, 1] x S), where S is the unit sphere of vectors 
u, 

2. evenness: g(—s, —uj) = g(s, to), 

3. moment conditions: for any integer k > 0, the fcth moment 



extends from the unit sphere S to a homogeneous polynomial of degree k in lo. 

The seemingly "trivial" evenness condition is sometimes the hardest to generalize to other 
transforms of Radon type, while it is often easier to find analogs of the moment conditions. 
This is exactly what happens in TAT. 

For the operators W, Ai in TAT, some sets of range conditions of the moment type had 
been discovered over the years [7,84,85,111], but complete range descriptions started to 
emerge only since 2006 [2, 4-6, 9, 43, 78]. 

Range descriptions for the more general operator W are harder to obtain than for Ai, 
and complete range descriptions are not known for even dimensions or for the case of the 
variable speed of sound. 

Let us address the case of the spherical mean operator Ai first. 




oo 




— oo 
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3.5.1 The range of the spherical mean operator M.. 

The support and smoothness conditions are not hard to come up with, at least when S is 
a sphere. By choosing appropriate length scale, we can assume that the sphere is of radius 1 
and centered at the origin, and that the interior domain Q is the unit ball B — {x \ \x\ — 1}. 
If / is smooth and supported inside B (i.e. / G Cq°(B)), then it is clear that the measured 
data satisfies the following 

Smoothness and support conditions: 

geC™(Sx [0,2]). (12) 

An analog of the moment conditions for g(y, r) := A4 f was implicitly present in [7, 84, 85] 
and explicitly formulated as such in [111]: 

Moment conditions: for any integer k > 0, the moment 

oo 

M k (y) = J r 2k+d - l g(y,r)dr (13) 
o 

extends from S to an (in general, non-homogeneous) polynomial Qk{x) of degree at most 2k. 

These two types of conditions happen to be incomplete, i.e. infinitely many others exist. 
The Radon transform experience suggests to look for an analog of evenness conditions. And 
indeed, a set of conditions called orthogonality conditions was found in [5,9,43]. 

Orthogonality conditions: Let — X 2 be the eigenvalue of the Laplace operator A in B with 
zero Dirichlet conditions and ipk be the corresponding eigenf unctions. Then the following 
orthogonality condition is satisfied: 

J g(x,t)d u Mx)jn/2-i(\t)t n - l dxdt = Q. (14) 

5x[0,2] 

Here j p {z) = c p z~ p J p {z) is the so called spherical Bessel function. 

The range descriptions obtained in 2D in [9] and then in general dimension in [5] showed 
that these three types of conditions completely describe the range of the operator M. on 
functions / G C^°(B). At the same time, the results of [5,43] showed that the moment 
conditions can be dropped in odd dimensions. It was then discovered in [2] that the moment 
conditions can be dropped altogether in any dimension, since they follow from the other two 
types of conditions: 

Theorem 25. [2] Let S be the unit sphere. A function g(y, t) on the cylinder S x IR + can be 
represented as M.f for some f G C^°(B) if an only if it satisfied the above smoothness and 
support and orthogonality conditions (12), (14). 

The statement also holds in the finite smoothness case, if one replaces the requirements 
byfe H* Q {B) and g G H s +{n ~ 1)/2 (S x [0,2]). 

The range of the forward operator M. has not been described when S is not a sphere, 
but, say, a convex smooth closed surface. The moment and orthogonality conditions hold 
for any S, and appropriate smoothness and support conditions can also been formulated, at 
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least in the convex case. However, it has not been proven that they provide the complete 
range description. 

It is quite possible that for non-spherical S the moment conditions might have to be 
included into the range description. 

A different range description of the Fredholm alternative type was developed in [103] (see 
also [45] for description of this result). 



3.5.2 The range of the forward operator W. 

We recall that the operator W (see (4)) transforms the initial value / in (3) into the 
observed on S values g of the solution. There exist Kirchhoff-Poisson formulas representing 
the solution p, and thus g = Wf in terms of the spherical means of / (i.e., in terms 
of M.f). However, translating the result of Theorem 25 into the language of W is not 
straightforward, since in even dimensions these formulas are non-local [30, 36] (pp. 682 and 
801 correspondingly): 

W(v,f) = ;s#-7ST TaT t n - 2 (Mf)(y,t), for odd n. (15) 



2V{n/2) \tdt 
and 

Wf(s , t) = _L_ ( "- 2) ' 2 [ r~HMf)(y,r) for _ 

Jyyi 1 T(n/2)\tdtJ J y/W^T 1 V ' 



The non-locality of the transformation for even dimensions reflects the absence of Huy- 
gens' principle (i.e. absence of sharp rear fronts of waves) in these dimensions; it also causes 
difficulties in establishing the complete range descriptions. In particular, due to the inte- 
gration in (16) A4f(y,t) does not vanish for large times t anymore. One can try to use 
other known operators intertwining the two problems (see [5] and references therein), some 
of which do preserve vanishing for large values of t, but this so far has lead only to very 
clumsy range descriptions. 

However, for odd dimensions, the range description of W can be obtained. In order to 
do so, given the TAT data g(y, t), let us introduce an auxiliary time- reversed problem in the 
cylinder B x [0, 2]: 

'qtt - Aq = for (x,t) G B x [0,2]), 
q(x,2) = q t (x,2) = for x E B, (17) 
q(y,t) =g(y,t) for (y, t) E S x [0, 2]). 

We can now formulate the range description from [43,45]: 

Theorem 26. [43, 45] For odd dimensions n and S being the unit sphere, a function g E 
C£°(S x [0,2]) can be represented as Wf for some f E C^°(B) if and only if the following 
condition is satisfied: 

The solution q of (17) satisfies qt(x, 0) = for all x E B. 

Orthogonality type and Fredholm alternative type range conditions, equivalent to the 
one in the theorem above, are also provided in [43,45]. 
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3.6 Speed of sound reconstruction 

Unsurprisingly, all inversion procedures outlined in Section 4 rely upon the knowledge of 
the speed of sound c(x). Although often, e.g. in breast imaging, the medium is assumed 
to be acoustically homogeneous, this is not a good assumption in many other cases. It has 
been observed (e.g., [68,71]) that replacing even slightly varying speed of sound with its 
average value might significantly distort the image; not only the numerical values, but also 
the shapes of interfaces between the tissues will be reconstructed incorrectly. Thus, the 
question of estimating c(x) correctly becomes important. One possible approach [71] is to 
use an additional transmission ultrasound scan to reconstruct the speed beforehand. The 
question arises of whether one could determine the speed of sound c(x) and the tomogram 
f(x) (assuming that / is not zero) simultaneously from the TAT data. In fact, one needs 
only to determine c(x) (without knowing /), since then inversion procedures of Section 4 
would apply to recover /. 

At the first glance, this seems to be an overly ambitious project. Indeed, if we denote 
the forward operator W by W c , to indicate its dependence on the speed of sound c(x), then 
the problem becomes, given the data g, to find both c and / from the equality 

Wcf = g- (18) 

A similar situation arises in the SPECT emission tomography (see [76, 88, 89] and references 
therein), where the role of the speed of sound is played by the unknown attenuation. It is 
known, however, that in SPECT the attenuation can be recovered for a "generic" /. 

What is the reason for such a strange situation? It looks like for any c one could solve 
the equation (18) for an /, and thus no information about c is contained in the data g. 
This argument is incorrect for the following reason: the range of the forward operator, as 
we know already from the previous section, has infinite co-dimension. Thus, this range has 
a lot of space to "rotate" when c changes. Imagine for an instance that the rotation is so 
powerful that for different values of c the ranges have only zero (the origin) in common. 
Then, knowing g in the range, one would know which c it came from. Thus, the problem of 
recovering the speed of sound from the TAT data is closely related to the range descriptions. 

Numerical inversions using algebraic iterative techniques (e.g., [143,144]) show that re- 
covering both c and / might be indeed possible. 

Unfortunately, very little is known at the moment concerning this problem. Direct usage 
of range conditions attempted in [68] has lead only to extremely weak and not practically 
useful results so far. A revealing relation to the transmission eigenvalue problem well known 
in inverse problems (see [29] for the survey) was recently discovered by D. Finch. Unfor- 
tunately, the transmission eigenvalue problem remains still unresolved. However, one can 
derive from this relation the following (still not too useful for TAT) uniqueness of the speed 
of sound determination result, due to M. Agranovsky: 

Theorem 27. If two speeds satisfy the inequality C\(x) > c 2 (x) for all x G Q and produce 
for some functions fi, f 2 the same non-zero TAT data g (i.e., VVci/i = g, Wc 2 / 2 = g), then 
ci(x) = c 2 (x). 

It is known [70, Corollary 8.2.3] that if a function f(x) is such that Af(x) ^ and for 
two acoustic speeds C\{x) and c 2 (x) it produces the same TAT data g, then c\ — c 2 . 
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It is clear that the problem of finding the speed of sound from the TAT data still requires 
significant analysis. 

4 Reconstruction formulas and procedures 

Numerous formulas, algorithms and procedures for reconstruction of images from TAT mea- 
surements have been developed by now. Most of these techniques require the data being 
collected on a closed surface (closed curve in 2D) surrounding the object to be imaged. Such 
methods are discussed in Section 4.1. We review methods that work under the assumption 
of constant speed of sound in Section 4.1.1. The techniques applicable in the case of the 
known variable speed of sound are considered in Section 4.1.2. Closed surface measurements 
cannot always be implemented, since in some practical situations the object cannot be com- 
pletely surrounded by the detectors. In this case, one has to resort to various approximate 
reconstruction techniques as discussed in Section 4.2. 

4.1 Full data (closed acquisition surfaces) 
4.1.1 Constant speed of sound 

When the speed of sound within the tissues is a known constant, the TAT problem can 
be reformulated (see Section 2) in terms of the values of the spherical means of the initial 
condition f(x). These means can be easily recovered from the measurements of the acous- 
tic pressure using formulas (15) and (16) (see the discussion in [7]). In this case, image 
reconstruction becomes equivalent to inverting the spherical mean transform Ai. Thus, in 
what follows, we consider the problem of reconstructing a function f(x) supported within 
the region bounded by a closed surface S from known values of its spherical integrals g{y,r) 
with centers on S: 

g(y,r)= j f{y + ru)r n - l duj, y e S, (19) 

where duo is the standard measure on the unit sphere. 

Series solutions for spherical geometry 

The first inversion procedures for the case of closed acquisition surfaces were described 
in [94,95], where solutions were found for the cases of circular (in 2D) and spherical (in 
3-D) surfaces, respectively. These solutions were obtained by the harmonic decomposition 
of the measured data and of the sought function f(x), followed by equating coefficients of 
the corresponding Fourier series. In particular, the 2D algorithm of [94] pertains to the case 
when the detectors are located on a circle of radius R. This method is based on the Fourier 
decomposition of / and g in angular variables: 

oo 

fix) = J2fk(p)e lkv , x = (pcos(^),psin(^)) (20) 

— oo 
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g{y(6),r) = ]>> fc (r)e^, y = (Rcos(0), Rsm(0)), 



— oo 



where 



/ 

Jo 



■oo 



(HmU) (s) = 2tt 



u(t)J m (st)tdt, 



is the Hankel transform and J m (t) is the Bessel function. As shown in [94], the Fourier 
coefficients fk{p) can be recovered from the known coefficients gk{v) by the following formula: 



This method requires division of the Hankel transform of the measured data by the 
Bessel functions Jk, which have infinitely many zeros. Theoretically, there is no problem: 
the range conditions (Section 3.5) on the exact data g imply that the Hankel transform 
H.q [(2irr)~ 1 gi t (r)} has zeros that cancel those in the denominator. However, since the mea- 
sured data always contain errors, the exact cancelation does not happen, and one needs a 
sophisticated regularization scheme to guarantee that the error remains bounded. 

This difficulty can be avoided (see, e.g. [78]) by replacing the Bessel function J in the 
inner Hankel transform by the Hankel function Hq. This yields the following formula for 



Unlike J m , Hankel functions Hm (t) do not have zeros for any real values of t, which removes 
the problems with division by zeros [94] . (A different way of avoiding divisions by zero was 
found in [62]) 

This derivation can be repeated in 3D, with the exponentials e lke replaced by the spherical 
harmonics, and with cylindrical Bessel functions replaced by their spherical counterparts. By 
doing this, one arrives at the Fourier series method of [95] (see also [135]). The use of the 
Hankel function Hq 1 ^ above is similar to the way the spherical Hankel function hjp is utilized 
in [95] to avoid the divisions by zero. 

Eigenfunction expansions for a general geometry 

The series methods described in the previous section rely on the separation of variables 
that occurs only in spherical geometry. A different approach was proposed in [82]. It 
works for arbitrary closed surfaces, but is practical only for those with explicitly known 
eigenvalues and eigenf unctions of the Dirichlet Laplacian in the interior. Such surfaces 
include, in particular, spheres, half-sp heres, cylinders, cubes and parallelepipeds, as well as 
the surfaces of crystallographic domains. 

Let and u m (x) be the eigenvalues and an ortho- normal basis of eigenfunctions of the 




fM ■ 
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Dirichlet Laplacian —A in the interior Q of a closed surface S: 



Au m {x) + \ 2 rn u m (x) = 0, 160, OCR", (21) 
u m (x) = 0, x e S, 

\\ u m\\l = J \u m (x)\ 2 dx = 1. 
h 

As before, one would like to reconstruct a compactly supported function f(x) from the 
known values of its spherical integrals g(y, r) (see (19)) with centers on S. Since u m {x) is the 
solution of the Dirichlet problem for the Helmholtz equation with zero boundary conditions 
and the wave number X m , this function admits the Helmholtz representation 

f d 

ujx) = J $ Am (|x - y\)— u m (y)ds(y) x e tt, (22) 

where §\ m (\x — y\) is a free-space Green's function of the Helmholtz equation (21), and n is 
the exterior normal to S. 

The function f(x) can be expanded into the series 

oc 

f( x ) = ^2 a mUm(x), where (23) 

m=0 

a m = u m (x)f(x)dx. 
Jq 

A reconstruction formula for a m (and thus for f(x)) will result, if one substitutes represen- 
tation (22) into (23) and interchanges the orders of integration: 

a m = j u m (x)f(x)dx = j I(y,X m )-^-u m (y)dA(x), (24) 

where 

I(y,\)= $ x (\x-y\)f(x)dx= g(y,r)$ x (r)dr. (25) 

Jn Jo 

Now f(x) can be obtained by summing the series (23). This method becomes computation- 
ally efficient when the eigenvalues and eigenf unctions are known explicitly, especially if a 
fast summation formula for the series (23) is available. This is the case for a cubic acqui- 
sition surface S, when the eigenf unctions are products of sine functions. The resulting 3D 
reconstruction algorithm is extremely fast and precise (see [82]). 

The above method has an interesting property. If the support of the source f(x) extends 
outside n, the algorithm still yields theoretically exact reconstruction of f(x) inside O. 
Indeed, the value of the expression (22) for all x lying outside Q is zero. Thus, when one 
computes (24) for x G M. n \ Q, values of f(x) are multiplied by zero and do not affect 
further computation in any way. This feature is shared by the time reversal method (see 
the corresponding paragraph in Section 4.1.2). The closed form FBP type reconstruction 
techniques considered in the next sub-section, do not have this property. In other words, 
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in presence of a source outside the measurement surface, reconstruction within f2 will be 
incorrect. 

The reason for this difference is that all currently known closed form FBP-type formulas 
rely (implicitly or explicitly) on the assumption that the wave propagates outside S in the 
whole free space and has no sources outside. On the other hand, the eig enfunction expansion 
method and the time reversal rely only upon the time decay of the wave inside S, which is 
not influenced by / having a part outside S. 



Closed form inversion formulas 

Closed-form inversion formulas play a special role in tomography. They bring about 
better theoretical understanding of the problem and frequently serve as starting points for 
the development of efficient reconstruction algorithms. A well known example of the use of 
explicit inversion formulas is the so-called filtered backprojection (FBP) algorithm in X-ray 
tomography, which is derived from one of the inversion formulas for the classical Radon 
transform (see, for example [73,88]). 

The very existence of closed form inversion formulas for TAT had been in doubt, till 
the first such formulas were obtained in odd dimensions by Finch et al in [42], under the 
assumption that the acquisition surface S is a sphere. Suppose that the function f(x) 
is supported within a ball of radius R and that the detectors are located on the surface 
S = dB of this ball. Then some of the formulas obtained in [42] read as follows: 
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/(*) 
/(*) 



8ir 2 R J \y-x\ V h 

dB 



8ir 2 R J \r dr 2 

OB 



1 d 2 



g(y,r] 



dA{y), 



r=\y—x\ 



1 d ( d g(y,r) 



8ir 2 R J \r dr \ dr r 

dB 



dA{y), 



(26) 
(27) 
(28) 



r=\y—x\ 



where dA(y) is the surface measure on dB and g represents the values of the spherical 
integrals (19). 

These formulas have a FBP (filtered back-projection) nature. Indeed, differentiation with 
respect to r in (27) and (28) and the Laplace operator in (26) represent the filtration, while 
the (weighted) integrals correspond to the backprojection, i.e. integration over the set of 
spheres passing through the point of interest x and centered on S. 

The so-called "universal backprojection formula" in 3D was found in [136] (it is also 
valid for the cylindrical and plane acquisition surfaces, see Section 4.2). In our notation, this 
formula takes the form 
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. 1 d g(y,r) 
n(y) -- 
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r dr 



dA{y), 



(29) 
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or, equivalently, 
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1 f d (I d g(y,r) 



dn \r dr r 



dB 



dA{y), 



(30) 



r=\y—x\ 



where n(y) is the exterior normal vector to dB. One can show [4, 90, 136] that formulas (26) 
through (29) are not equivalent on non-perfect data: the result will differ if these formulas 
are applied to a function that does not belong to the range of the spherical mean transform 
A4. A family of inversion formulas valid in M n for arbitrary n > 2 was found in [81]: 

f( x ) = ,/o_\n-i div J n (y) h (y^ \ x - y\) dA (v), ( 31 ) 

dB 



4(2tt) 



where 



h(y,t)= / Y(Xt) 



2R 



2R 



J(Xr)g(y,r)dr — J(Xt) / Y(Xr)g(y,r)dr 



X 2n ~ 3 dX, 



J(t) 



J, 



n/2-l 



(t) 



Y(t) 



Y n/2 -l{t) 



(32) 



(33) 



fn/2-l ' v 1 pi/2-X ' 

and J n /2-i{t) and Y n / 2 -i{t) are respectively the Bessel and Neumann functions of order 
n/2 — 1. In 3D, J(t) and Y(t) are simply t~ l sint and t~ x cost and formulas (31) and (32) 
reduce to (30). 

In 2D, equation (32) also can be simplified [4], which results in the formula 

_ 2R 
dB Lo 

where dB now stands for the circle of radius R and dl(y) is the standard arc length. 

A different set of closed-form inversion formulas applicable in even dimensions was found 
in [41]. Formula (34) can be compared to the following inversion formulas from [41]: 

2R 

f( x ) = 7T7^ A / / 9(y, r ) lo g0 2 -\x- y\ 2 ) dr dl(y), (35) 



dl(y), 



(34) 
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r— — ) log(r 2 — \x — y\ 2 ) dr dl{y). 



2txR J J dr \ dr r 

dB o 



(36) 



Finally, a unified family of inversion formulas was derived in [90]. In our notation, it has 
the following form: 
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K n (y,t) 



X 2n - 3 Y(Xt) I / J(Xr)g(y,r)dr I dX 
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<y - x,y - £> 
\x - y\ 



dA{y), 



(37) 
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where dB is the surface of a ball in W 1 of radius R, functions J and Y are as in (33), and £ 
is an arbitrary fixed vector. In particular, in 3D 



J(t) 



2 sint 



,J(t) 



2 cost 



71 t V 7T t 

and, after simple calculation, the above inversion formula reduces to 
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dA{y). 



(38) 



Different choices of vector £ in the above formula result in different inversion formulas. For 
example, if £ is set to zero, the ratio <y "j"^ / "j"^ > equals Rcosa, where a is the angle between 
the exterior normal n(y) and the vector y — x; when combined with the derivative in t this 
factor produces the normal derivative, and the inversion formula (38) reduces to (30). On 
the other hand, the choice of £ = x in (38) leads to a formula 
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r=\x—y\ 



which is reminiscent of formulas (26)-(28). 



Greens' formula approach and some symmetry considerations 

Let us suppose for a moment that the acoustic detectors could measure not only the 
pressure p(y, t) at each point of the acquisition surface S, but also the normal derivative 
dp/dn on S. Then the problem of reconstructing the initial pressure f(x) becomes rather 
simple. Indeed, one can use the knowledge of the free-space Green's function for the wave 
equation and invoke the Green's theorem to represent the solution p(x, t) of (3) in the form 
of integrals over 5" involving p(x, t) and its normal derivative and the Green's function and 
its normal derivative. (This can be done in the Fourier or time domains.) This would require 
infinite observation time, but in 3D the time T(Q) will suffice, afte r which the wave escapes 
the region of interest (a cut-off also would work approximately in 2D. similarly to the time- 
reversal method). This Green's function approach happens to be, explicitly or implicitly, 
the starting point of all closed form inversions described above. The trick is to rewrite the 
formula in such a way that the unknown in reality normal derivative dp/dn disappears from 
the formula. 

This was achieved in [81] by reducing the question to some integrals involving special 
functions and making the key observation that the integral 

I x (x,y) = j J(X\x-z\)^-Y(X\y-z\)dA(z), x,yeBcR n 

dB 

is a symmetric function of its arguments: 

I x (x, y) = I x {y, x) for x, y e B C R. n (39) 
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Similarly, the derivation of (37) in [90] employs the symmetry of the integral 
K x (x,y) = [ J(X\x-z\)Y(X\y-z\)dA(z), x,y e B cR n . 



dB 



In fact, the symmetry holds for any integral 



W x (x,y) = / U(X\x-z\)V(X\y-z\)dA(z), 



x,y G B c R n , 



dB 



where i7(A|x|) and V(A|x|) are any two radial solutions of Helmholtz equation 



Au(x) + X 2 u(x) = 0. 



(40) 



It is straightforward to verify this symmetry when S is a sphere and B is the corresponding 
ball, and the points x, y lie on the boundary S only, rather than anywhere in B. This follows 
immediately from the rotational symmetry of S. The same i s true for the normal derivatives 
on S of W\(x, y) in x and y. 

This boundary symmetry happens to imply the needed full symmetry (39) for x,y G B. 

Indeed, W\(x, y) is a solution of the Helmholtz equation separately as a function of x and 
of y. Let us introduce a family of solutions {w n (x)}^ =0 of (40) in B, such that the members 
of this family form an orthonormal basis for all solutions of the latter equation in B. For 
example, the spherical waves, i.e. the products of spherical harmonics and Bessel functions, 
can serve as such a basis. 

Then W\(x,y) can be expanded n the following series: 



Since W\(x,y) is a solution to the Helmholtz equation in dB x dB, coefficients b n ^ m are 
completely determined by the boundary values of W\. Since the boundary values are sym- 
metric, the coefficients are symmetric, i.e. b n)Tn = b m n which by (41) immediately implies 
W\(x, y) = W\(y, x) for all pairs (x,y) G B x B . 

This consideration extends to infinite cylinders and planes. This explains why the "uni- 
versal backprojection formula" (30) is valid also for infinite cylinders and planes [136]. Since 
the sort of symmetry used is shared only by these three surfaces, we believe it is unlikely 
that a closed-form formula could exist for any other acquisition surface. 

Algebraic iterative algorithms 

Iterative algebraic techniques are among the favorite tomographic methods of reconstruc- 
tion and have been used in CT for quite a while [73, 88, 89]. They amount to discretizing the 
equation relating the measured data with the unknown source, followed by iterative solu- 
tion of the resulting linear system. Iterative algebraic reconstruction algorithms frequently 
produce better images than those obtained by other methods. However, they are notori- 
ously slow. In TAT, they have been used successfully for reconstructions with partial data 
([16,17,108]), see Section 4.2. 
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Parametrix approaches 

Some of the earlier non-iterative reconstruction techniques [75] were of approximate na- 
ture. For example, by approximating the integration spheres by their tangent planes at 
the point of reconstruction and by applying one of the known inversion formulas for the 
classical Radon transform, one can reconstruct an approximation to the image. Due to the 
evenness symmetry in the classical Radon projections (see Section 3.5), the normals to the 
integration planes need only fill a half of a unit sphere, in order to make possible the recon- 
struction from an open measurement surface. A more sophisticated approach is represented 
by the so-called "straightening" methods [115, 116] based on the approximate reconstruction 
of the classical Radon projections from the values of the spherical mean transform Aif of 
the function f(x) in question. These methods yield not a true inversion, but rather what 
is called in micro-local analysis a parametrix. Application of a parametrix reproduces the 
function / with an additional, smoother term. In other words, the locations (and often the 
sizes) of jumps across sharp material interfaces, as well as the whole wave front set WF(f), 
are reconstructed correctly, while the accuracy of the lower spatial frequencies cannot be 
guaranteed. (Sometimes, the reconstructed function has a more general form Af, where A 
is an elliptic pseudo-differential operator [66, 121, 126] of order zero. In this case, the sizes of 
the jumps across the interfaces might be altered.) Unlike the approximations resulting from 
the discretization of the exact inversion formulas (in the situations when such formulas are 
known), the parametrix approximations do not converge, when the discretization of the data 
is refined and the noise is eliminated. Parametrix reconstructions can be either accepted 
as approximate images, or used as starting points for iterative algorithms. See [123] for a 
recent discussion of parametrices. 

These methods are closely related to the general scheme proposed in [22, 32] for the 
inversion of the generalized Radon transform with integration over curved manifolds. It 
reduces the problem to a Fredholm integral equation of the second kind, which is well suited 
for numerical solution. Such an approach amounts to using a parametrix method as an 
efficient pre- conditioner for an iterative solver; the convergence of such iterations is much 
faster than that of algebraic iterative methods. 

Numerical implementation and computational examples. 

By discretizing exact formulas presented above, one can easily develop accurate and effi- 
cient reconstruction algorithms. The 3D case is especially simple: computation of derivatives 
in the formulas (26)- (30) and (38) can be easily done, for instance by using finite differences; 
it is followed by the backprojection (described by the integral over dB), which requires pre- 
scribing quadrature weights for quadrature nodes that coincide with the positions of the 
detectors. The backprojection step is stable; the differentiation is a mildly unstable opera- 
tion. The sensitivity to noise in measurements across the formulas presented above seems 
to be roughly the same. It is very similar to that of the widely used FBP algorithm of 
classical X-ray tomography [88,89]. In 2D, the implementation is just a little bit harder: 
the filtration step in formulas (34)- (36) can be reduced to computing two Hilbert transforms 
(see [78]), which, in turn, can be easily done in the frequency domain. 

The number of floating point operations (flops) required by such algorithms is determined 
by the slower backprojection step. In 3D, if the number of detectors is m? and the size of 
the reconstruction grid is m x m x m, the backprojection step (and the whole algorithm) 
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will require 0(m 5 ) flops. In practical terms this amounts to several hours of computations 
on a single processor computer for a grid of size 129 x 129 x 129. 

In 2D, the operation count is just 0(m 3 ). As it is discussed in Section 2.4, the 2D 
problem needs to be solved, when integrating line detectors are used. In this situation, the 
2D problem needs to be solved m times in order t o reconstruct the image, which raises the 
total operation count to 0(m 4 ) flops. 

Figure 10 shows three examples of simulated reconstruction using formula (34). The 
phantom we use (Figure 10(a)) is a linear combination of several characteristic functions of 
disks and ellipses. Part (b) illustrates the image reconstruction within the unit circle from 257 
equi-spaced projections each containing 129 spherical integrals. The detectors were placed 
on the concentric circle of radius 1.05. The image shown in Figure 10(c) corresponds to the 
reconstruction from the simulated noisy data that were obtained by adding to projections 
values of a random variable scaled so that the L 2 intensity of the noise was 15% of the 
intensity of the signal. Finally, Figure 10(d) shows how application of a smoothing filter (in 
the frequency domain) suppresses the noise; it also somewhat blurs the edges in the image. 




Figure 10: Example of a reconstruction using formula (34): (a) phantom; (b) reconstruction 
from accurate data; (c) reconstruction from the data contaminated with 15% noise; (d) 
reconstruction from the noisy data with additional smoothing 

4.1.2 Variable speed of sound 

The reconstruction formulas and algorithms described in the previous section work under the 
assumption that the speed of sound within the region of interest is constant (or at least close 
to a constant). This assumption, however, is not always realistic - for example, if the region 
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of interest contains both soft tissues and bones, the speed of sound will vary significantly. 
Experiments with numerical and physical phantoms show [68, 71] that if acoustic inhomo- 
geneities are not taken into account, the reconstructed image might be severely distorted. 
Not only the numerical values could be reconstructed incorrectly, but so would the material 
interface locations and discontinuity magnitudes. 

Below we review some of the reconstruction methods that work in acoustically inhomo- 
geneous media. We will assume that the speed of sound c(x) is known, smooth, positive, 
constant for large x, and non-trapping. In practice, a transmission ultrasound scan can be 
used to reconstruct c(x) prior to thermoacoustic reconstruction, as it is done in [71]. 

Time reversal 

Let us assume temporarily that the speed of sound c is constant and the spatial dimension 
is odd. Then Huygens' principle guarantees that the sound wave will leave the region of 
interest Q in time T = cj (diamfi), so that p(x, t) = for all x e fl and t > T. Now one can 
solve the wave equation back in time from t = T to t = in the domain Q x [T, 0] , with zero 
initial conditions at T and boundary conditions on S provided by the data g collected by 
the detectors. Then the value of the solution at t — will coincide with the initial condition 
f(x) that one seeks to reconstruct. Such a solution of the wave equation is easily obtained 
numerically by finite difference techniques [52,68]. The required number of floating point 
operations is actually lower than that of methods based on discretized inversion formulas 
((9(m 4 ) for time reversal on a grid m x m x m in 3-D versus 0(m 5 ) for inversion formulas), 
which makes this method quite competitive even in the case of constant speed of sound. 

Most importantly, however, the method is also applicable if the speed of sound c(x) is 
variable and/or the spatial dimension is even. In these cases, the Huygens' principle does not 
hold, and thus the solution to the direct problem will not vanish within dfl in finite time. 
However, the solution inside Q will decay with time. Under the non-trapping condition, 
as it is shown in (11) (see [34,129,130]), the time decay is exponential in odd dimensions, 
but only algebraic in even- dimensions. Although, in order to obtain theoretically exact 
reconstruction, one would have to start the time reversal at T = oo, numerical experiments 
(e.g., [68]) and theoretical estimates [67] show that in practice it is sufficient to start at the 
values of T when the signal becomes small enough, and to approximate the unknown value 
of p(x, T) by zero (a more sophisticated cut-off is used in [123] , which leads to an equation 
with a contraction operator). This works [52,68] even in 2D (where decay is the slowest) 
and in inhomogeneous media. However, when trapping occurs, the "invisible" parts blur 
away (see Section 3.3 for the discussion). 

Eigenfunction expansions. 

An "inversion formula" that reconstructs the initial value f{x) of the solution of the 
wave equation from values on the measuring surface S can be easily obtained using time 
reversal and Duhamel's principle [3]. Consider in Q the operator A = — c 2 (x)A with zero 
Dirichlet conditions on the boundary S = dQ. This operator is self-adjoint, if considered in 
the weighted space L 2 (Q; c~ 2 (x)). Let us denote by E the operator of harmonic extension, 
which transforms a function (ft on S to a harmonic function on Q which coincides with (ft on 
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S. Then / can be reconstructed [3] from the data g in (3) by the following formula: 



oo 

f{x) = iEg\ t=0 ) - j A~hmirA-2)Eig tt )ix,r)dT, (42) 
o 

which is valid under the non-trapping condition on c(x). However, due to the involvement 
of functions of the operator A, it is not clear how useful this formula can be. 

One natural way to try to implement numerically the formula (42) is to use the eigen- 
function expansion of the operator A in fl (assuming that such expansion is known). This 
quickly leads to the following procedure [3]. The function fix) can be reconstructed inside 
Q from the data g in (3), as the following L 2 (l?)-convergent series: 

fix) = J2fkM%), (43) 
k 

where the Fourier coefficients f k can be recovered from the data using one of the following 
formulas: 

00 

fk = K 2 9k(0) - Xf f sin (\ k t)g»(t)dt, 
o 

oo 

fk = K 2 9kiO) + A^ 2 J cos (X k t)g' k (t)dt, or (44) 
o 



/* = ~ X k 1 / sin (\ k t)g k (t)dt = -A^ 1 / / sin (X k t)g(x, t)^-(x)dxdt, 

OS 



where 



g(x,t)^-(x)dx. 



One notices that this is a generalization of the expansion method of [82] discussed in 
Section 4.1.1 to the case of a variable speed of sound. Unlike the algorithm of [82], the 
present method does not require the knowledge of the whole space Green's function for A 
(which is in this case unknown). However, computation of a large set of eigenfunctions and 
eigenvalues followed by the summation of the series (43) at the nodes of the computational 
grid may prove to be too time consuming. 

It is worthwhile to mention again that the non-trapping condition is crucial for the 
stability of any TAT reconstruction method in acoustically inhomogeneous media. As it was 
discussed in Section 3.4, trapping can significantly reduce the quality of reconstruction. It 
is, however, most probable that trapping does not occur much in biological objects. 



4.2 Partial (incomplete) data 

Reconstruction formulas and algorithms of the previous sections work under the assumption 
that the acoustic signal is measured by detectors covering a closed surface S that surrounds 
completely the object of interest. However, in many practical applications of TAT, detectors 
can be placed only on a certain part of the surrounding surface. Such is the case, for 
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example, when TAT is used for breast screening - one of the most promising applications of 
this modality. Thus, one needs methods and algorithms capable of accurate reconstruction 
of images from partial (incomplete) data, i.e. from the measurements made on open surfaces 
(or open curves in 2D). 

Most exact inversion formulas and methods discussed above are based (explicitly or 
implicitly) on some sort of the Green's formula, Helmholtz representation, or eigenfunction 
decomposition for closed surfaces, and thus they cannot be extended to the case of partial 
data. The methods that do work in this situation rely on approximation techniques, as 
discussed below. 

4.2.1 Constant speed of sound 

Even the case of an acoustically homogeneous medium is quite challenging when reconstruc- 
tion needs to be done from partial data (i.e., when the acquisition surface S is not closed). 
As it was discussed in Section 3.3, if the detectors located around the object in such a way 
that the "visibility" condition is not satisfied, accurate reconstruction is impossible: the 
"invisible" interfaces will be smoothed out in the reconstructed image. On the other hand, 
if the visibility condition is satisfied, the reconstruction is only mildly unstable (similarly to 
the inversion of the classic Radon transform) [103,123]. If, in addition, the uniqueness of 
reconstruction from partial data is guaranteed (which is usually the case, see Section 3.3.1), 
one can hope to be able to develop an algorithm that would reconstruct quality images. 

Special cases of open acquisition surfaces are a plane or an infinite cylinder, for which 
exact inversion formulas are known (see, for example, [18,40,48,92,138] for the plane and 
[139] or for a cylinder). Of course, the plane or a cylinder would have to be truncated in any 
practical measurements. The resulting acquisition geometry will not satisfy the visibility 
condition, and material interfaces whose normals do not intersect the acquisition surface will 
be blurred. 

Iterative algebraic techniques (see the corresponding paragraph in Section 4.1.1) were 
among the first methods successfully used for reconstruction from surfaces only partially 
surrounding the object (e.g., [16, 17, 108]). As it is mentioned in Section 4.1.1, such methods 
are very slow. For example, reconstructions in [17] required the use of a cluster of computers 
and took 100 iterations to converge. 

Parametrix type reconstructions in the partial data case were proposed in [19]. A couple 
of different parametrix-type algorithms were proposed in [105,107]. They are based on 
applying one of the exact inversion formulas for full circular acquisition to the available 
partial data, with zero-filled missing data and some correction factors. Namely, since the 
missing data is replaced by zeros, each line passing through a node of the reconstruction grid 
will be tangent either to one or to two circles of integration. Therefore some directions during 
the backprojection step will be represented twice, and some only once. This, in turn, will 
cause some interfaces to appear twice stronger then they should be. The use of weight factors 
was proposed in [105, 107] in order to partially compensate for this distortion. In particular, 
in [105] smooth weight factors (depending on a reconstruction point) are assigned to each 
detector in such a way that the total weight for each direction is exactly one. This method is 
not exact; the error is described by a certain smoothing operator. However, the singularities 
(or jumps) in the image will be reconstructed correctly. As shown by numerical examples in 
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[105], such a correction visually significantly improves the reconstruction. Moreover, iterative 
refinement is proposed in [105, 107] to further improve the image, and it is shown to work 
well in numerical experiments. 

Returning to non-iterative techniques, one should mention an interesting attempt made 
in [111, 112]) to generate the missing data using the moment range conditions for A4 (see 
Section 3.5). The resulting algorithm, however, does not seem to recover the values well; 
although, as expected, it reconstructs all visible singularities. 

An accurate 2D non-iterative algorithm for reconstruction from data measured on an 
open curve S was proposed in [83]. It is based on pre-computing approximations of plane 
waves in the region of interest Q by the single layer potentials of the form 

Z(X\y -x\)p{y)dl(y), 



where p{y) is the density of the potential, which needs to be chosen appropriately, dl(y) 
is the standard arc length, and Z(t) is either the Bessel function Jo(t), or the Neumann 
function Y Q (t). Namely, for a fixed £ one finds numerically the densities p^ t j{y) and p^y{y) 
of the potentials 



Wj(x,p u ) = / J (\\y - x\)pt,j(y)dl(y), (45) 
Js 

W Y (x,ps, Y )= f Y (\\y-x\)pt, Y (y)dl(y), (46) 
Js 

where A = |£|, such that 

Wj(x, p^j) + Wy(x, p^y) ~ ex P( - ^ ' x ) f° r an x efl. (47) 

Obtaining such approximations is not trivial. One can show that exact equality in (47) cannot 
be achieved, due to different behavior at infinity of the plane wave and the approximating 
single-layer potentials. However, as shown by numerical examples in [83] , if each point in Q is 
"visible" from S, very accurate approximations can be obtained, while keeping the densities 
p^j and p^y under certain control. 

Once the densities p^j and p^ Y have been found for all £, function f(x) can be easily 
reconstructed. Indeed, for the Fourier transform /(£) of f(x) 

/(0 = TT" / f( x ) ex P(-^ • x )dx, 



2tt _ 

one obtains, using (47) 

m)*Y I f(x) [Wj(x,p t j)+W Y (x,p t Y)}dx 



l 

1 



S 



f(x)J (X\y - x\)dx 
f(x)Y (\\y - x\)dx 



Pz,j(y) dl (y) 

Pt, Y (y)dl(y), (4* 



39 



where the inner integrals are computed from the data g: 



f(x)J (X\y - x\)dx = / g(y,r)J (Xr)dr, 

J R+ 

f{x)Y (X\y - x\)dx = / g(y,r)Y (Xr)dr. 
n Jr+ 



(49) 
(50) 



Formula (48), in combination with (49) and (50), yields values of /(£) for arbitrary £. 
Now f(x) can be recovered by numerically inverting the Fourier transform, or by a reduction 
to a FBP inversion [73, 88] of the regular Radon transform. 

The most computationally expensive part of the algorithm, which is computing the den- 
sities p^j and p^Y, needs to be done only once for a given acquisition surface. Thus, for a 
scanner with a fixed S, the resulting densities can be pre-computed once and for all. The 
actual reconstruction part then becomes extremely fast. 

Examples of reconstructions from incomplete data using this technique of [83]) are shown 
in Figure 11. The images were reconstructed within the unit square [—1, 1] x [—1,1], while 
the detectors were placed on the part of the concentric circle of radius 1.3 lying to the left 
of line :ri = l. We used the same phantom as in Figure 10(a)); the reconstruction from the 
data with added 15% noise is shown in Figure 11(b); part (c) demonstrates the results of 
applying additional smoothing filter to reduce the effects of noise in the data. 





Figure 11: Examples of reconstruction from incomplete data using the technique of [83]. 
Detectors are located on the part of circular arc of radius 1.3 lying left of the line X\ — 1. (a) 
reconstruction from accurate data (b) reconstruction from the data with added 15% noise 
(c) reconstruction from noisy data with additional smoothing filter 



4.2.2 Variable speed of sound 

The problem of numerical reconstruction in TAT from the data measured on open surfaces 
in the presence of a known variable speed of sound currently remains largely open. One 
of the difficulties was discussed in Section 3.3: even if the speed of sound c(x) is non- 
trapping, it can happen that some of the characteristics escape from the region of interest to 
infinity without intersecting the open measuring surface. Then stable reconstruction of the 
corresponding interfaces will become impossible. It should be possible, however, to develop 
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stable reconstruction algorithms in the case when the whole object of interest is located in 
the visible zone. 

The generalization of the method of [83] to the case of variable speed of sound is so 
far problematic, since this algorithm is based on the knowledge of the open space Green's 
function for the Helmholtz equation. In the case of a non-constant c(x), this Green's func- 
tion is position-depended, and its numerical computation is likely to be prohibitively time- 
consuming. 

A promising approach to this problem, currently under development, is to use time 
reversal with the missing data replaced by zeros, or maybe by a more clever extension (e.g., 
using the range conditions, as in [111, 112]). This would produce an initial approximation to 
f(x), which one can try to refine by fixed-point iterations; however, the pertinent questions 
concerning such an algorithm remain open. 

An interesting technique of using a reverberant cavity enclosing the target to compensate 
for the missing data is described in [31]. 

5 Final remarks and open problems 

We list here some unresolved issues of mathematics of TAT/PAT, as well as some develop- 
ments that were not addressed in the main text. 

1. The issue of uniqueness acquisition sets S (i.e., such that transducers distributed along 
S provide sufficient information for TAT reconstruction) can be considered to be re- 
solved, for most practical purposes. However, there remain significant unresolved the- 
oretical questions. One of them consists of proving an analog of Theorem 8 for non- 
compactly supported functions with a sufficiently fast (e.g., super-exponential) decay 
at infinity. The original (and only known) proof of this theorem uses microlocal tech- 
niques [7, 122] that significantly rely upon the compactness of support. However, one 
hopes that the condition of a fast decay should suffice for this result. In particular, 
there is no proven analog of Theorem 6 for non-closed sets S (unless S is an open part 
of a closed analytic surface). 

Techniques developed in [42] (see also [8] for their further use in TAT) might provide 
the right approach. 

This also relates to still unresolved situation in dimensions 3 and higher. Namely, one 
would like to prove Conjecture 9. 

2. Concerning the inversion methods, one notices that closed form formulas are known 
only for spherical, cylindrical, and planar acquisition surfaces. The question arises 
whether closed form inversion formulas could be found for any other closed surface? It 
is the belief of the authors that the answer to this question is negative. 

Another feature of the known closed form formulas that was mentioned before is that 
they do not work correctly if the support of the sought function f(x) lies partially 
outside the acquisition surface. Time reversal and eigenfunction expansion methods 
do not suffer from this deficiency. The question arises whether one could find closed 
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form formulas that reconstruct the function inside S correctly, in spite of it having 
part of its support outside. Again, the authors believe that the answer is negative. 

3. Besides algebraic iterative approaches, there are no reliable reconstruction methods in 
the case of the detectors partially surrounding the target, if the medium is acoustically 
inhomogeneous (see Section 4.2.2). This contrasts with the acoustically homogeneous 
situation (Section 4.2.1). 

4. The complete range description of the forward operator W in even dimensions is still 
not known. It is also not clear whether one can obtain complete range descriptions 
for non-spherical observation sets S or for a variable sound speed. The moment and 
orthogonality conditions do hold in the case of a constant speed and arbitrary closed 
surface, but they do not provide a complete description of the range. For acoustically 
inhomogeneous media, an analog of orthogonality conditions exists, but it also does 
not describe the range completely. 

5. The problem of unique determination of the speed of sound from TAT data is largely 
open. 

6. As it was explained in the text, knowing full Cauchy data of the pressure p (i.e., its 
value and the value of its the normal derivative) on the observation surface S leads to 
unique determination and simple reconstruction of /. However, the normal derivative 
is not measured by transducers and thus needs to be either found mathematically or 
measured in a different experiment. Thus, feasibility of techniques [13,28] relying on 
full Cauchy data requires further mathematical and experimental study. 

7. In the standard X-ray CT, as well as in SPECT, the so called local tomography 
technique [37-39, 77] is often very useful. It allows one to emphasize in a stable way 
singularities (e.g., tissue interfaces) of the reconstruction, even in the case of incomplete 
data (in the latter case, the invisible parts will be lost). An analog of local tomography 
can be easily implemented in TAT, for instance, by introducing an additional high- 
frequency filter in the FBP type formulas. 

8. The mathematical analysis of TAT presented in the text did not take into account 
the issue of modeling and compensating for the acoustic attenuation. This subject 
is addressed in [24,74,87,113,120], but probably cannot be considered completely 
resolved. 

9. The initial pressure f(x) that was the center of all discussions in the chapter (as well 
as in most papers devoted to TAT/PAT), is related, but is not exactly identical to 
the optical features of interest of the tissue. The issue of recovering the actual optical 
parameters of the tissue after the initial pressure f(x) is found is non-trivial and is 
addressed, probably for the first time, in [20]. 

10. This chapter as well as most other papers devoted to TAT/PAT is centered on the 
initial pressure f{x). This quantity is related, but is not exactly identical to the 
relevant optical features of the tissue. The problem of recovering the actual optical 
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parameters of tissue (after f(x) is found) is non-trivial and is addressed, probably for 
the first time, in [20]. 

11. The TAT technique discussed in the chapter uses active interrogation of the medium. 
There is a discussion in the literature of a passive version of TAT, where no irradiation 
of the target is involved [110]. 
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